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■ Abstract. Hybrid Monte Carlo simulations that implement the fermion action using multiple terms 

^ I are commonly used. By the nature of their formulation they involve multiple integration time scales 

. in the evolution of the system through simulation time. These different scales are usually dealt with 

by the Sexton-Weingarten nested leapfrog integrator. In this scheme the choice of time scales is 
somewhat restricted as each time step must be an exact multiple of the next smallest scale in the 
sequence. A novel generalisation of the nested leapfrog integrator is introduced which allows for far 
greater flexibility in the choice of time scales, as each scale now must only be an exact multiple of 
, the smallest step size. 
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INTRODUCTION 



Hybrid Monte Cario(HMC)[l] is the algorithm of choice for generating lattice gauge 
\^ '. field configurations that include the effects of fermion loops. The main expense in 

O I such simulations is evaluating the contribution from the fermion determinant. A variety 

of improvements to the basic HMC algorithm have been developed in an effort to 
ameliorate this expense. Many of these variants involve the introduction of additional 
terms into the fictitious Hamiltonian. For example, Clark and Kennedy[2] introduce 
multiple pseudofermion fields into the Rational HMC algorithm, each contributing a 
fraction of the fermion determinant. Alternatively, Hasenbuch[3] uses the fermion matrix 
rS \ (with a heavier quark mass) as a preconditioner to the desired fermion action. This causes 

I the fermion action to be split into two parts. Another example is the polynomial filtered 

HMC algorithm[4], which uses short polynomial approximations to the inverse fermion 
matrix in order to separate the dynamics of the Hamiltonian. Multiple levels of filtering 
can be introduced, and this technique is also applicable to single flavour simulations, 
allowing the number of terms to be simulated in the Hamiltonian to grow to five or 
more. The Schwarz preconditioning technique employed by Luscher[5] also makes use 
of a hierarchy of time scales. 

The traditional way of dealing with a multiple time scale integration is to use the 
nested leapfrog algorithm by Sexton and Weingarten[6]. This scheme is quite restrictive 
however, as beginning with the finest time scale, the step- size for the next (coarser) time 
scale must be an exact multiple of the previous scale. This may prevent one from being 
able to choose the most efficient set of integration parameters, particularly if there are 
many time scales. Here we present an improvement of this scheme, which allows for 



greater freedom of choice, as each time scale must only be a multiple of the finest scale. 



Hybrid Monte Carlo Overview 

In this section we provide a brief overview of the Hybrid Monte Carlo algorithm[l] 
in order to provide a framework for the introduction of our novel integrator. We wish 
to generate an ensemble {Ui} of representative gauge fields distributed according to 
the probability distribution p{Ui) = e^'^t^'l, where the effective action for full QCD 
S[U] = Sg[U] + Sf[U] can be divided into two parts, the gauge action Sg[U] and the 
fermion action Sf[U]. Here we have assumed that the fermionic degrees of freedom 
have been integrated out in the usual way. 

In the Hybrid Monte Carlo algorithm, the quantum lattice field theory is embedded 
in a higher-dimensional classical system through the introduction of a fictitious (simu- 
lation) time. The gauge field U is associated with its (fictitious) conjugate momenta P, 
and the classical system is described by the Hamiltonian, 

H[U,P] = T[P]+S[U], (1) 

where the fictitious kinetic energy is given by T[P] = Y.x,^ jT^^Pfi{x)^- Given a config- 
uration U, a new gauge field U' is generated by performing a HMC update U — j- U', 
which consists of two steps: 

(1) Molecular Dynamics Trajectory: Sample P from a Gaussian ensemble. Integrate 
Hamilton's equations of motion to deterministically evolve {U ,P) along a phase 
space trajectory to {U',P'). 

(2) Metropolis step: Accept or reject the new configuration {U',P') with probability 

p{U U') = mm{l,e-^),AH = H[U',P'] -H[U,P]. 

The discretised equations of motion are derived by requiring that the Hamiltonian be 
conserved along the phase space trajectory. We can express the equations of motion in 
terms of the time evolution operators induced by the kinetic and potential energy terms. 
The evolution operators that evolve the gauge field and its conjugate momenta forward 
a simulation time step h are given respectively by 

Vrih) : {U,P} ^ {Ucxp {ihP),P}, (2) 

Vs{h):{U,P}^{U,P-hU^}. (3) 

Note that we demand that the evolution must preserve the SU (3) property of the gauge 
field. 

We must then combine these evolution operators into an overall evolution operator 
that is reversible. The simplest such integration scheme is the leapfrog 

yH{h)=Vs{^)VT{h)Vs{\). (4) 

After discretisation, for sufficiently small step sizes h, the integration will conserve 
the Hamiltonian up to 0{h^). 



Multiple time-scales in molecular dynamics integrators 



If the action and thus the Hamikonian is split into two parts Hi and H2, 

H = T[P]+^Si[U] +S2[U] (5) 

Hi H2 

then we define integrators for Hi and H2 as follows 

VHi{h)=Vs,{i)VT{h)Vs,{l), VH,{h)=Vs,{h). (6) 

A compound integrator for the full Hamiltonian can be constructed by using a Sexton- 
Weingarten scheme[6]: 



VH{h)=VH, 



■T 



m 



^//.(^) (7) 



where m G Z. This nested leapfrog integrator effectively introduces two time-scales into 
the evolution, h and hjm. Additional time scales may be introduced by repeating the 
nesting procedure. 

RESULTS 

We begin by assuming that our Hamiltonian consists of at least three terms, 

// = r + 5i+52 + ..., (8) 

where T is the "kinetic energy" due to the conjugate momenta, and the terms Si imple- 
ment the lattice QCD action. Typically we would choose Si = 5g, that is is the gauge 
action. Then 52,53, .. . will be the terms implementing the fermion action according to 
our algorithm of choice, be it mass-preconditioning or polynomial filtering and so on. 
The only thing we assume about the fermion action is that it is implemented such that the 
only fields that are updated during the integration are the gauge field U and its conjugate 
momenta P (e.g. using pseudofermions). 

For each scale we associate a timestep hi and a corresponding integer A^, such that 
hi = l/Ni. We require that i = 1 corresponds to the scale at which the gauge field is 
updated. 

Assuming that A^, > Nj for i < j a nested leapfrog algorithm then requires that 

A^,|A^,_iV/>l. (9) 

That is, A^,- must be a divisor of A'^,- 1 , or equivalently, hi must be an exact multiple of 1 . 
This means, for example, that each successive scale must be at least twice as coarse as 
the previous scale. It also may lead to being forced to choose a scale that is smaller 
than the one desired for a given time scale in order to simultaneously control the finite 
step- size errors as well as maintain the required arithmetic relation (9). The restriction 
of having to choose successive divisors for the various Ni may not be the most efficient 
or flexible way of performing the molecular dynamics integration. 



A generalised leap-frog integrator 



We present a generalised integration scheme in which it is only required that the 
integration scales satisfy the relation 

Ni\Ni\/i>l. (10) 

This is of course equivalent to requiring that the step size hi is an exact multiple of hi. 

In a standard leapfrog algorithm, one alternates between updates Vt to the gauge field 
U and updates V5 to the conjugate momenta. Let V,- denote the update to P corresponding 
to the action 5, . Now, as the guide bosons are held fixed during an integration the updates 
Vi only depend upon the gauge field. As the updates Vj are additive to P, it follows that 
the different V,- commute: 

v,{^)Vj{hjM'^) = Vj{hj)Vi{hi) = Vi{hi)Vj{hj). (11) 

Define the integers 

mi = Ni^Ni (12) 

to be the ratios of the scales. In order to construct our reversible integrator we first define 
a map 

, (V if m\k 

@[V;m,keN] = < .[ . , ^ . (13) 
[ /(the identity) otherwise. 

Let niT be the lowest common multiple of {m,}, and let hj be the smallest time step (in 
our case hx = hi). Then our integrator is 



k=l 



where h = inxhT is the total timestep taken by V. The above expression is straightfor- 
wardly implemented in code. We demonstrate this with a pseudocode implementation 
here. Denote hy {a = b mod m} the usual notion of congruence modulo m. Then we 
can implement the generalised integrator as follows. 

Pseudo-code for the generalised integrator: 

• For each term in the action 5,- perform an initial half-step Vi{^hi) updating P. 

• Loop over j = I to N — I 

- Apply VV(/?) to update U. 

- If {0 = j mod m,} apply V'/(/?/) to update P 

• Apply Vrih) to update U. 

• For each term in the action 5, perform a final half-step Vi{hhi) updating P. 



The advantage of the generalised integrator is that it allows finer control over the 
different scales. An analysis of the finite-step size errors for the generalised integrator is 
provided in the next section. 



Integrator Error Analysis 

We perfomi an error analysis of our generalised integrator for a simple choice of step- 
sizes, following the procedure in [6]. Given a Hamiltonian H we can write the evolution 
operator for our system as V{h) = exp (hH), with step-size h. Here we have defined H 
as the linear operator on the vector space of functions / on phase space {p,q) defined 
by the Poisson bracket 

y\dpidqi dqjdpij' 

If we write the Hamiltonian as 

H = T + Si+S2 + S3 + S4 + ... (16) 

then for each term in the Hamiltonian we can correspondingly define a linear operator 
using the Poisson bracket relation above. We denote the linear operator associated with 
a given term in the Hamiltonian by adding a "hat" to the appropriate symbol. 

Proceeding with the error analysis, we make use of the Baker-Campbell-Hausdorff 
result. 



.XA^XB^XA ^ (^^2A + B) + ^([[A,5],A] + [[A,B],B]) + 0{X')) (17) 

and apply this to the generalised leapfrog integrator in the simple case of H = T + Si + 
52, where the time scale for each term in H corresponds to a number of integration steps 
Nt — 6jNi = 3 and N2 = 2 respectively. The integrator for this simplest non-trivial case 
can be written as 



, ^ ha ha h 'T ho h 'T ha h 'T ho h 'T ha ha 

Vnih) = e^-e^^'e^^e^^'e^^e^^^e^^e^^'e^^e^^'e^^^. (18) 

Repeated application of our BCH result allows us to deduce that the above expression 
can be written as 



Vh (h) = exp (/iH + (1 [[52, Tlf] + ^ [[52, f] , 52] + 



:^[[SuT],Si] + —[[SuT]J] + -[[Si,T],S2])) (19) 

From this expression we can immediately see that the error in the generalised integrator 
relative to the leading term is C>(/z^), just as for the regular leapfrog. 
If we examine the individual leapfrog integrators corresponding to 

Hi^T + Su H2^T + S2 (20) 

we obtain 

, ^ /7 c h T' ha h T' ha h 't ha 

VHiih) = e6^ie3^e3^ie3^e3^ie3^e6^i (21) 
= exp {hHi +h\^[[Si,f],f] + J^[[5i, f ],5i])), (22) 



and 



Vn^ih) = e4^2e5^e5^2e5^eT^2 (23) 

= exp (^hH2 + h\^[[S2M + ^[[^2,r],52])) (24) 

Hence we see that the only difference between the individual integrators and our 
generalised integrator is the cross term [[5i,7'],52]. The algorithm is identical to a 
standard nested leapfrog in the case where Ni\Ni-i. 

CONCLUSIONS 

We have introduced a novel integration scheme that generalises the nested leapfrog 
scheme by Sexton and Weingarten[6]. The new scheme has the advantage that each 
term in the fictitious Hamiltonian may be assigned a time step that only need be an exact 
multiple of the finest time step. This is an improvement on the nested leapfrog, where 
each scale must be an exact multiple of the next smallest scale in the hierarchy. 

Each term in the fictitious Hamiltonian will have a corresponding "force term" as- 
sociated with it. The typical size of this force term leads one to choose an appropriate 
integration time scale for that term. The large the force term, the smaller the time scale 
that is required to keep the finite step- size errors under control. With the novel integra- 
tion scheme introduced here one has more freedom to choose a time-scale that is appro- 
priate to the force associated with a given term - in the nested leapfrog one may have 
been forced to choose a time scale that was smaller than needed due to the restrictions 
imposed by the size of the next smallest scale in the hierarchy. 

In a HMC simulation with two degenerate flavours where the fermion action has been 
split into two (or more) pieces, we already have at least three different time scales: one 
for the gauge action, and two (or more) for the fermion action. With the development of 
multiple time scale algorithms such as polynomial filtered HMC [4] that are applicable to 
both two flavour and single flavour simulations, the number of terms in the Hamiltonian 
can grow quite large. It is in these cases that the novel scheme proposed here will become 
particularly useful. 
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